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Abstract 

A new algorithm to calculate Coulomb wave functions with all of its arguments complex is 
proposed. For that purpose, standard methods such as continued fractions and power/ asymptotic 
series are combined with direct integrations of the Schrodinger equation in order to provide very 
stable calculations, even for large values of \r]\ or Moreover, a simple analytic continu- 

ation for ^{z) < is introduced, so that this zone of the complex z-plane does not pose any 
problem. This code is particularly well suited for low-energy calculations and the calculation of 
resonances with extremely small widths. Numerical instabilities appear, however, when both 
\ri\ and are large and |3fJ(^)| comparable or smaller than 

Program Summary 

Title of the program: cwfcomplex 
Catalogue number: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 
Program summary URL: 
Licensing provisions: none 

Computers on which the program has been tested: DELL GX400 

Operating systems: Linux, Windows 

Programming language used: C++ 

Memory required to execute with typical data: 

No. of bits in a word: 64 



^E-mail address: nmichel@utk.edu 
Phone : 1-865-576-4295 



Fax : 1-865-576-8746 



1 



No. of processors used: 1 



Has the code been vectorized?: No 

No. of bytes in distributed program, including test data, etc.: 
No. of lines in distributed program: 2422 
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1 Introduction 

Coulomb wave functions are one of the most basic objects of particle theory. They describe the 
behavior of a particle in a point-like Coulomb field, and thus appear in virtually all domains of 
quantum physics. The correspondent dimensionless Coulomb Schrodinger equation reads: 



where w{z) is a Coulomb wave function, ^ is the orbital angular momentum, and rj the Som- 
merfeld parameter. 

The Coulomb wave functions can be expressed with hyper-geometric functions [1]. The 
regular Coulomb wave function reads: 






(2) 



Q(r7) = 2' exp 



TIT] + [log(r(l + £ + If])) + log(r(l + £ - IT]))] 
2 



log(r(2£ + 2)) (3) 
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In this expression, uo can be equal to ±1 and the normahzing Gamow factor C^(?7) [I] is given. 
Outgoing {uo = 1) and incoming {uo = —1) Coulomb wave functions are defined the following 
way: 

H^^^{z) = e-[--''i°g(2-)-'^§+-^('')] ,Fo (^-i + tuor], l + i + luor]; ; (4) 

iog(r(i + £ + 277))-iog(r(i + f-zr/)) 
Mv) = 2i 

where the Coulomb phase shift (Te{r]) appears \T\. The analytic continuation for complex i 
and rj of Ref. [21 [3] for the function log(r(z)) occurring in Ce{ri) and cre{r]) is followed, thus 
guaranteeing consistent values even when the negative real axis branch cut of complex variables 
1 + £ + if] and 1 + £ — ir] is crossed. The regular Coulomb wave function F^^^, as well as the 
logarithmic irregular Coulomb wave function Ge^rj, can be expressed with if^^ and if^"^ [1]: 

Despite the deceptively simple form of Eq. ([T]) and analytical expressions of Coulomb wave 
functions of Eqs. (|2|^ . the Coulomb wave function is difficult to compute numerically. Already 
on the real axis, it can vary by many orders of magnitude for moderate values of \ri\. The 
situation becomes even worse when the wave function is analytically continued to the complex 
plane. Analytic continuation arises when one deals, for example, with resonant states, as 
energies become complex [4j. It appears also with non-integer values of £ with, for example, 
Regge pole trajectory calculations [3]. Coulomb wave functions are multivalued functions of the 
complex variable z in the general case and thus a branch cut must be imposed on the negative 
z-real axis [2J. This implies that numerical calculations must be employed with care, as wave 
functions issued directly from standard numerical expressions do not follow the same branch 
cut discontinuities as the requested Coulomb wave function. 

The Coulomb wave function computation has been considered in many papers. A recent 
review of numerical methods and definitions for both non-relativistic and relativistic cases can 
be found in Ref. [B]. Most of papers have dealt with only real arguments [7|, or with at least 
one of them real {£ in Ref. [8], z in Ref. [9]). The special important case of Whittaker functions, 
with rj purely imaginary (bound Coulomb wave functions), has also been treated on its own [TOj. 
The first paper (and only one to our knowledge) which considered all complex arguments in a 
unified way is Ref. pij. Through the use of continued fractions calculated with the powerful 
Lentz method, as well as recurrence relations in £ and Pade approximants, the authors managed 
to encompass a large part of the complex plane for each £, r], and z. The program of Ref. [11] 
quickly became a standard in the physics community and is part of the CERNLIB library [12] . 
However, important parts of the complex plane remained uncovered for both numerical and 
theoretical reasons, as was already stated in [TT]. First, because of numerical instabilities of 
used recursions, one cannot calculate Coulomb wave functions, for example, close to imaginary 
axes when the modulus of rj or Q{£) becomes large. As a consequence, the case of large \r]\ and 
integer £ for all z, important for low-energy narrow resonant states, has remained problematic 
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jl3j . Moreover, due to the different branch cuts of Coulomb wave functions and asymptotic 
series/continued fractions in the complex ^-plane, it is impossible to directly calculate H^^^{z) 
for ^{z) < and uQ{z) < 0. 

In order to circumvent these caveats, it has been chosen to complement standard methods 
with direct integrations of Eq. ([1]). The latter can be simply implemented, and the only require- 
ment is that one has to integrate in directions of increasing modulus of the wave function to 
avoid numerical instability Also, the ^{z) < and ujQ{z) < parts of the complex 2;-plane 



can be accessed with numerical methods, as the (wrong) Coulomb wave functions coming out 
of the latter are linear combinations of the true Coulomb wave functions whose coefficients can 
be computed precisely, so that their determination becomes straightforward. With these new 
features, it will be demonstrated that the range of arguments is much larger than in previous 
programs. It is also important to state that quadruple precision is not needed with the proposed 
method. 

The structure of the paper is as follows: first, the used numerical methods will be described 
in Sec. ([2]). Examples of calculations will then be depicted for several sets of arguments in 
Sec. (|3]). In particular, the determination of resonant states with extremely small widths will 
be discussed. The structure of the program will then be described. Finally, conclusions and 
perspectives will be stated. 

2 Numerical methods 
2.1 Power series for 

The regular solution Fi^^j can be expanded in power series pQ: 

+ 00 

Fe,,{z) = Civ) ^a„;^"+^+^ (8) 

n=0 

ao = 1 

V 



i+1 

2r] a^-i - an-2 
n{n + 2i+l) 



\fn>2 (9) 



This formula is very useful for small values of \z\, but is unstable for large \z\ because of 
numerical cancellations. Hence, it is used only for \z\ < -. 



4 



2.2 Asymptotic series 



i?^^ can be expanded in asymptotic series, so that for \z\ large enough, H^^^{z) can be calculated 
up to a given numerical precision with a finite number of terms [15] : 



N-l , 

Or. 



Hl^{z) - e'-[-^iog(2^)-^f+-.W] ^ ^ (10) 

n=0 ^" 

On+i = TTl — rT\ ^" Vn > 11 

2iuj[n + 1) 

However, due to the different branch cuts discontinuity of if^^ and the asymptotic series, 
Eq. flTU]) is correct only for Iki^z) > or 3fJ(z) < and uj^{z) > ^llj. In the rest of the 
complex plane, it nevertheless provides a linear combination of Ht and H7, which is utilized 
to determine if^^ (see Sec. (12. 6p ). 

In practice, the asymptotic series give a meaningful result for a given if |aAr/2;^| < e with 



e the numerical precision [15]. In addition, one checks if the Wronskian of the two functions 
generated by Eq. (fTOl) respectively using u = 1 and = — 1 is equal to 2i. The Wronskian 
value can be evaluated with Eq. ( JTOl) for \z\ —>■ +oo. If 3^(2;) > 0, Fi^r}{z) is calculated with 
asymptotic series using Eq. ([6]) if H/^^{z) and H^,^{z) are correctly computed with Eq. ffTOl) . 



2.3 Continued fractions 



The logarithmic derivatives = F^^/Fi^r^ and h"^ = can be expanded in continued 

fractions [H]: 



h-{z) 



•+i . 1 

\-tUJ + - 

z z 



1-^ +- 



—liujaz —2iuj{a + l)z 

b + 2iuoz+ h+l + 2iuoz ^ 

ac (a + l)(c + 1) 

2{z -r] + iuj)+ 2{z-ri + 2iuj) + 



(12) 
(13) 



where the standard notations a = 1 + £ + iurj, b = 2i + 2, and c = —i + iur] are used pjj. The 
value of (also denoted as /) is derived from Eq. ([2]) and is thus theoretically independent of 
ui. Lentz method is used to evaluate continued fractions numerically [TTj . 

The domain of convergence is the whole complex plane besides F^^^ zeros, while the one 
of h'^ follows 2-^0 analytic properties, so that it is the whole complex plane minus the half-axis 
[0 : —iuJOo[, where h'^ has a branch cut discontinuity. 

The continued fraction h"^ is particularly important, as with the knowledge of Fe^ri{z), F^^{z), 

W: 



and the Wronskian relation F^'^if^^ 



1, it can be used to determine H'^^^{z) 



FeM [f(z)-h-(z)] 



(14) 



Note that H'^^{z) and F^r]{z) must be numerically linearly independent for this formula to be 
stable. If they are not, H^!^{z) is calculated instead and H'^^^{z) can be deduced from it and 
Eq. ([6]). As h'^ and if^^ have different branch cuts, h'^ is equal to the logarithmic derivative of 
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H^^ only if '^{z) > or '^{z) < and uj'^{z) > 0, so that Eq. ( fT4l) is correct in this zone only. 
However, the continued fraction can be used even outside this zone if one takes care of branch 
cuts (see Sec. (12. 6p ). Added to that, the continued fractions and h'^ play a prominent role 
for the calculation of Coulomb wave functions by direct integration (see Sec. (I2.4p ). 

The numerical applicability of these continued fractions is, however, hindered by spurious 
effects. It has been noticed in Ref. [16] that Eq. f|T2|) exhibits anomalous convergence. When 
1^1 becomes large, the general term of /'^(z) becomes very small before increasing very much, 
and then only to decrease again to provide a convergent result. As a consequence, both values 
of f~^{z) and f~{z) are always calculated and compared to check convergence. However, the 
anomalous convergence phenomenon is weaker when one chooses u such that Lv'is{z) < [16], 
so that f^{z) can be correct even if f^{z) ^ f {z) numerically. 

The case of is much better [Hj, but problems have nevertheless been encountered. For 
example, the numerical value of h^{z) for £ = 0, ?7 = 10, and z = 0.01 — 3i is wrong and 
numerically equal to f{z). This difficulty is removed by using only the uj for which \ f{z) — h'^{z)\ 
is large enough (i.e. larger than 1 or at least larger than \ f{z) — h^'^{z)\). 

Another problem is the very slow convergence of Eq. f|T3l) in the vicinity of the branch cut 
for moderate \z\ (see Table ([T])). Direct integration is used to solve this problem. For that, if 
the number of iterations in Lentz method exceeds 100,000 (one also assumes 5R(z) > 0), one 
calculates H'^^^^zq) and H'^^^{zo) with zq not too close to the imaginary axis, chosen so that 
increases from zo to z. The slow convergence of h'^ is absent for zq, so that Eq. (fT3|) can be used 
for the integration starting point. Then, one integrates Eq. ([1]) from zq to z, which is a stable 
operation as increases along the integration path, h'^i^z) is then equal to H^^^i^z) / H'^^^{z) 

at the end of integration. If 3fJ(z) < 0, one uses the symmetry formula h';^^^{z) = —hj^'^^{—z). 
This formula can be demonstrated using the fact that H'^^^{z) oc Hl^^^{—z) for uQ{z) > (both 
functions are solutions of Eq. ([1]) and are minimal in the considered region for \z\ +oo), and 
analytic continuation as both functions fi:z^ h'^,rj{^) fi'-z^ —hj'^^^—z) have the same 
branch cut. 



2.4 Direct integration 

Considering the simplicity of the Coulomb equation (Eq. ([1])), direct integration is a suitable 
method to calculate Coulomb wave functions. For that, the Burlisch-Stoer-Henrici method 
of Ref. [13] is used. However, one has to pay attention to two problems. Firstly, no branch 
cut discontinuity can come out of direct integration, so that it is necessary to integrate in the 
zones of the complex plane where branch cut effects are absent. Hence, numerical integration 
is performed only for ^{z) > 0. For the other half of the complex plane, one uses the symmetry 
transformation z —z, rj —7], leaving Eq. ([1]) invariant. Secondly, numerical integration is 
stable only if the modulus of the Coulomb wave function increases or remains close to constant. 
Increase or decrease of the wave function along the integration path is determined by its second- 
order Taylor expansion at z = zq + h: 

where is either F^^^ or -ff^^, h the integration step, and {zq, "^{zq), "^'{zq)) the starting point of 
the numerical integration. If the modulus of the ratio defined in Eq. f|T5l) is larger than one, the 
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numerical integration can be performed safely. If not, the continued fraction q{z) = '$'{z)/'${z) 
is evaluated with Eq. (IT^ (g = f^) or Eq. (IT^ (g = h"^). Eq. ([T]) is integrated backward 
from z to Zq with {z, 1, q{z)) as the starting point, guaranteeing stable integration. One 
then obtains {zq, "^dzo), "^'^{zq)) after integration, with obviously ^'c(^o) = ^(^o)/^(^) and 
"^'^{zo) = '^'{zo)/'^{z). The value of {z, "^{z), "^'{z)) comes forward. The only nuisance in 
this method is that the continued fraction q{z) can be wrong due to numerical instability (see 
Sec. (E3D). 

This can be partially solved if one considers a direct integration of Fn^^^z). If 3ft(£) > —1, F^^^ 
increases in modulus with |2;| in the vicinity of z = (non-oscillatory zone). As a consequence, 
if \Fi^ri\ is found to decrease on its initial path, zq is reinitialized to z/{2\z\), where the power 
series formula of Eq. ([8]) is available, so that a decrease of |-Ff,r)| from zq to z is less likely to 
happen. If the direct integration of Fe^ri{z) is found to be unstable despite this change of path, 
it is preferred to calculate Hf^^{z) with direct integration {u chosen so the branch cut of /i^^ is 
avoided) as h'^ is numerically more stable than (see Sec. (12. 3p ). Fi^riiz) is then calculated 
with Eq. ([6]) and the following formula: 

^''^^""^ ^ Hl^{z) [h^{z)-h~-{z)] ^^^^ 

which can be obtained similarly to Eq. f|T^ . This process is, however, stable if H'^^{z) and 
H^^{z) are not numerically equal, so that it is not employed if l-F^,,,] is found to be smaller 
than 0.1 on its integration path. This is a sound procedure as f^lz) is usually correct in this 
case if u'^i^z) < (see Sec. (12. 3p ). 

As a result, direct integration is a powerful tool to determine Fi^ri{z) outside the zone of 
applicability of Eq. ([8]). It also provides h^{z) close to its branch cut, where the numerical cost 
of Eq. ( fT3l) becomes prohibitive (see Sec. (12.31) ). 

2.5 H^^ expansion 

When 2£ is not an integer, the following formula can be used to calculate if^^ [TT] : 



= ^^^^>- , —'-^'^ (17) 
' ' sm X 

X = a.iri) - ^-i-iiv) - + l/2)7r (18) 

In practice, it has been chosen to apply it only for |Q'(i')| > 1 and \z\ < 1, as other methods 
have been found to be more robust for other cases. For a given z, the expression of Eq. f[T7|) is 
numerically stable if the Wronskian relation between and F_i_i^^ is respected: 

Fl^F_e-i,r, - i^A.^Vi,, = sinx (19) 

Note that Eq. f|T8l) is not used to calculate sinx as it is unstable due to cancellation effects. 
Another formula is preferred: 

sinx = -(2£+ 1) Ce{7]) C^e-i{v) (20) 
which can be demonstrated using Eqs. flHfTOl) with 2; — > 0. 
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2.6 Analytic continuation for ^(z) < 

Analytic continuation for 3?(z) < is first considered for the regular function Ff^jj- Using 
Eqs. Fe^niz) and F^, -ri{—z) can be shown to be proportional: 

Fe,r,{z) = -e-^(^-''^ F,,_^{-z) for arg(^) > 

Fe,r,{z) = -e--(''+^^) F,^_^{-z) for arg(^) < (21) 

Hence, Fi^riiz) can always be deduced from z), so that calculations for > are 

sufficient to determine Fi^^] in all the complex plane. 

The situation is more complicated for -f^^^, as the direct evaluation of Eq. ffTOl) and Eq. f|T^ 
provides correct values for 3^^(^) < and uj'^{z) > 0, but wrong results occur when ^{z) < 
and uj'^{z) < We will denote as H'^J'^^''^ and H^^^^'^'^ the numerical values coming from 
a naive implementation of respectively Eq. ( fTOl) and Eq. (fT4l) . As they are issued from analytic 
expressions providing solutions of Eq. ([T]), they are still solutions of this equation even when 
< and uj^{z) < 0. However, the different branch cuts of Hl^, H';'J'^^''\ and H';'J^^'^ 
imply that the two latter functions are linear combinations of H^^ and H^^ in this quadrant of 
the 2;-complex plane. Their coefficients will be shown to be very simple expressions of i and rj 
and can be related to standard circuital relations [T7] . 

If one considers \z\ +oo, one can deduce from Eqs. fHHTOj) that: 

H'^^^iz) = H'^J^''\z) + ^^'''\z) for ^{z) < and uj%z) < (22) 

where a^^ is a constant depending on i and rj. The equality H^^ = H^^ ^^^"^^ in the considered 
region is also used. From Eqs. (I61I22I) . one has: 

2zFe,,{x.) = Hl^^'^\x^) + (a+ - 1) Hl^^^''^\x^) (23) 
-2^F,,,(x+) = Hl^^^''^\x^) + (a_ - 1) <(^^'^^(x+) (24) 

where x+ = x + ie and X- = x — ie for x < and e > 0. Branch cut discontinuities of F^^^ and 
H'lJ'^^'^^ are straightforward from Eqs. ( PfTOj) . so that Eq. ( l23l) can be rewritten as: 



2z F,,,(x+) = e-'-^ hI^^^''\x+) + e'^^ (a+ - 1)H;J'''''\x+) + 0(e) (25) 

Finally, using Eqs. fl24|25p with e — >• 0, one obtains a^^ = 1 — e^*'^^*''"^'^) and hence the requested 
formula: 

H^^iz) = Hi;^^'\z) + [1 - e2-(*'^-M] ^^^'\z) (26) 

for which '^{z) < and uj'^{z) < 0. 

Considering the different branch cut discontinuities of Eqs. fllO|14l) on the negative real axis 
and analytic continuation, one obtains: 

H;J^'^'\z) =e^'^^'''~'''^Hi;^^'\z) (27) 

with ^{z) < and uj'^{z) < 0. 
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Using Eqs. (16|26II27I) . the formulas analog to Eq. (126|) for continued fractions are derived: 



H'^^^iz) = H';f^'\z) - 2iuj [e^-^^-*") - l] (28) 
Hl^{z) = ^''^''\z) + 2tuj e-2-(^-+^^) (29) 

with K(z) < 0, ujQ{z) < for Eq. ([28]) but u'^iz) > for Eq. ([29]). As the calculation of 
Fg^riiz) is prerequisite to determine H^^{z) with continued fraction formulas (see Eq. f[T^ ). the 
numerical evaluation of H'fj^{z) with Eqs. ( !28|29l) is straightforward. 

Even though the expressions of the coefficients in front of H^^ in Eq. ( l26l) and Fi^^j in 
Eqs. ( I28f29l) are elementary, care must be given to calculate them due to possible overflow or 
underflow and numerical cancellations. One has to use complex generalizations of the standard 
C-language functions loglp{x) = log(l + x) and expml{x) = — 1 for x ^ to avoid possible 
numerical inaccuracies. 

As a result, the 3^(2;) < domain no longer poses any numerical problem, as the formulas 
of Eqs. (I26ll28f29l) render it comparable to the rest of the 2;-complex plane. 



2.7 Poles of the Coulomb wave functions 

When 1 + / + iuT] is a negative integer. Coulomb wave functions are undefined (see Eqs. ( 13]|5]) ). 
Nevertheless, if one considers 3?(£) > — 1, numerical solutions of Eq. ([1]) always exist and can be 
computed. For this, Fe^rj is defined with Eq. ([8]) putting arbitrarily Ci{ri) = 1. Direct integration 
can be performed precisely for F^^^, as the continued fraction of Eq. (|T2l) is finite so that no 
numerical inaccuracy can occur. H^^ can still be defined with Eq. IHM as h~'^ ^ f^, so that 
one can calculate two linearly independent solutions of Eq. ([1]) when 1 + / + iurj is a negative 
integer. Note that the branch cut of H^^ is, for this definition, [0 : iu!oo[ and not the negative 
real axis. H'^^ and Ge^n can, however, not be defined so that they are arbitrarily put equal to 
Fi^jj and H^^ respectively. 



2.8 Quasi- real £, j] and z 



When i, 1] or z are very close to their real axes with at least one of them complex, the imaginary 
part of Coulomb wave functions can become tens of order of magnitude smaller than their real 
parts. Consequently, it can be numerically imprecise as calculations are always provided up to 
the same absolute precision for both real and imaginary parts. This is especially visible if one 
deals with resonant states of extremely small widths 7 such as proton emitters (7 ~ 10~^° keV) 
[18] . For these kinds of states, one generally uses approximate current formulas [19], providing 
very good values for 7 if it is small enough. However, it is possible to reach the same precision 
directly at the Coulomb wave function level. For this, one expands the Coulomb wave functions 
Fi^^ and Gi^rj up to first order in the vicinity of the real axes of i, t], and z: 



A{i,ri,z) = A{ir,T]r,x) + i 



OA dA OA 



(30) 



where Iri^lriX and ii,r]i,y are respectively the real and imaginary parts of i, rj, and z, and 
A{i,r],z) is either Fe^riiz) or Ge^n{z). All values involving A{ir,VrjX) are real (one considers 
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X > only) so that function and partial derivatives can be evaluated numerically. In practice, if 
e is the demanded numerical precision, the conditions \y\ < y/emm{l,x), \rii\ < and \£i\ < y/e 
must be fulfilled for Eq. (150]) to be used. It was checked that direct and approximate results 
yield the same results up to a precision comparable to e if \y\ = y/emm{l,x), = \rii\ = -y/e 
with e = 10-1° (see Table ((2])). 

H'^j^{z) can be obtained straightforwardly from the knowledge of Fi^^{z) and G'^^^(z). Con- 
sequently, the whole Coulomb wave function can be derived up to a given relative numerical 
precision for both real and imaginary parts even if they are very different in modulus. 



2.9 Scaled wave functions and alternative normalization 

It often happens that Hf^{z) overflows or underflows when \z\ or \ri\ become large, but only 
through the exponential factor of Eq. (jl]). As a consequence, the following scaled Coulomb 
wave functions can also be calculated in the code: 

Hl^,, = HIM e-^-l^-^'-s^^^)] (31) 

They are particularly useful if one calculates products of Coulomb wave functions where the 
different exponential factors cancel each other. 

At the limit of very small energies, where |?7| is very large, Ci{rj) can also overflow or 
underflow, so that it is no longer possible to calculate Coulomb wave functions. However, their 
normalization factor is usually unimportant, as in the case of a resonant state calculation. For 
that, we introduced the following renormalized wave functions: 

F. ^M" = Ciivr' F,,M 
HU^y = Ceiv) HIM 

G, ,M" = Ciiv) Ge,r,iz) (32) 
An example of a resonant state for which Ci{ri) underflows will be given in Sec. (13.21) . 



2.10 Recurrence relations and associated Wronskian tests 

Coulomb wave functions obey recurrence relations of their angular momentum i [TT] : 

S 1 

where w is any of the F,G,H^ or functions. Re = (2i + 1) and Se = - + ^. 

3fJ(£), denoted ir, is supposed to be larger than zero. These recurrence relations are stable 
provided |w£,,,(2;)| increases with £r. As Fe^j^{z) —>■ with £r +oo [H], for ir large enough, 
the recurrence relations are stable with ir decreasing if one calculates regular Coulomb wave 
functions and with ir increasing if one calculates irregular Coulomb wave functions. For the 
irregular wave functions, one has the most stable calculations if one calculates H^^^{z) with uj 
chosen so \H'^^^^^{z)\ < |_fZ'^-'^(2;)|, io being the angular momentum of smallest modulus. Indeed, 



10 



this guarantees if^ ^ to be the minimal solution of Eq. ([T]) if Fi^^^i is not. One may have, 
however, a turning point It before which \F^^ri{,z)\ increases [H]. In this case, F^^r) has to be 
recurred backward from the angular momentum of largest modulus, denoted £i, backward to 
It but forward from £o to If Conversely, if^^ must be recurred backward from It to £o and 
forward from It to £i. 

The previous recurrence relations provide additional relations between Coulomb wave func- 
tions: 

Ftr^H^^^ - Fi^i^r^H"^ = — — (34) 

If a; is chosen so < \Hl!^{z)\, in order to have two Coulomb wave functions numerically 

linearly independent, Eqs. fl34ll35p provide a good test to check the accuracy of the Coulomb 
wave functions calculated with the methods of previous sections, as angular momentum recur- 
rence relations do not enter them. 

The non-standard normalization defined in Eq. (132!) can also be used in the code along with 
recurrence relations, which can be obtained straightforwardly from Eqs. (I32f33p . Scaling of H"^^ 
for u = ±1 (see Eq. (!3T!) )) is, however, not considered in this context, because all Coulomb 
wave functions have to be numerically finite for the method to work, so that it would only 
result in the trivial multiplication of H^^[z) by e~*'^[^~'''°s*^^^)l. 

3 Examples 

3.1 Calculations in difficult zones of the complex plane 

In order to illustrate the proposed numerical methods, we selected sets of £, t], and z parameters 
having sizable values, z is always of the form RtC^^ , where Rt = \ri\ + -|- 1)| -|- is a 

generalization of the turning point in the complex plane. Problems of convergence indeed 
typically occur in the vicinity of this point [11]. \r]\ and \i\ have been chosen so that at least 
one is large in a set of parameters (see Figs. (fT][2][3]|4l[5!) ). 

It was chosen not to have both \r]\ and large with |3?(^)| < |3(^)| or comparable, as 

calculations become unstable therein. For these values of £, ?7, and 2, Coulomb wave functions 
vary by several orders of magnitude along the complex circle of radius Rt and oscillate much as 
well. Hence, we represented the decimal logarithm of the modulus of Coulomb wave functions, 
which varies smoothly. One should note that the discontinuities encountered at ^ = vr follow 
the branch cuts imposed to Coulomb wave functions and are not induced by any numerical 
inaccuracy. These calculations show that Coulomb wave functions can be calculated precisely 
even when they vary much in argument and modulus. 

3.2 Calculations of resonant states of very small widths 

In order to show the possibilities of the present program related to resonant states with ex- 
tremely small widths, e.g. proton emitters we consider a spherical Schrodinger equation 
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with a Woods-Saxon potential crudely mimicking a heavy nuclear target acting on a proton 
projectile: 



h=^ ^° „ , + (36) 

2/i i+exp(^) ^ ^ ^ 

where /i is the reduced mass of the proton state so that h'^/2fi = 21 MeV fm^, d is the diffuseness 
of the potential fixed at 0.63 fm, Rq is the radius of the potential of 6.5 fm, Vq is the depth of 
the Woods-Saxon potential, and Vc{r) is a Coulomb potential generated by a uniformly charged 
sphere of radius Rq and charge Z = 66. These parameters correspond to the proton emitter 
i^^Ho [19]. 

The 2si/2 proton state energy and width of the Hamiltonian h of Eq. is calculated 
for several values of Vq (see Table ([3])). The width 7 coming from direct integration of h is 
compared with the following standard current approximation 



ic = m) - 
/i 



u{R) 



H+{kR) 



(37) 



where k is the linear momentum of the proton resonant state, u{r) its radial wave function, 
and R a radius large enough so u{r) oc H^{kr) for r > R. The value = 20 fm was chosen. 
As expected, both values 7 and 7c are identical, as 7 is very small. If Vq = 56.46 MeV, one 
obtains an energy of 6.949-10"^ MeV while the width is numerically zero. This value could be 
computed only through the renormalization of Coulomb wave functions of Eq. fl32|) . as for this 
energy log;^Q(C£(?7)) = —535, implying Cg^rj) underflow. 

These results show that the proposed program is very well suited for the direct calculation of 
very narrow resonances, for which one has to enter numerically challenging areas of the complex 
plane. It would be interesting to use this program along with coupled-channel integration 
methods [181 [19] in order to extend the current method to deformed states. 



4 The program cwf complex 
4.1 Routines of the program 

The code cwfcomplex is written in standard C++, uses only standard libraries and is thus 
portable on many machines. It is separated in four different files: complex-functions. H, cwf- 
comp.H, cwfcomp.cpp and tesLrec-rel.cpp. 

complex Junctions. H contains elementary complex functions which are not in the standard 
library, and routines calculating constants specific to the Coulomb wave functions: 

• inf_norm: provides the infinite norm of a complex number. 

• isfinite: returns true if the complex number is finite. 

• operators overloading of complex and integers. 

• expml : complex generalization of the function expml{x) = — 1 precise for x 0. 
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• loglp: complex generalization of the function loglp{x) = log(l + x) precise for x — > 0. 

• log_Gamma\ log(r(z)) calculated with the method of Ref. [3]. 

• sigmaJ-calc: complex Coulomb phase shift ceiv) (see Eq. ([5])). 

• log_CLeta_calc: log(C^(?7)) (see Eq. ([3])). 

• log_cut-ConstanLAS-calc: logarithm of the constant in front of H^!^ (-^^d) (126|) . 

• log_cut_constanLCFa_calc: logarithm of the constant in front of F^^^ in Eq. ( l28l) . 

• log_cut_constant_CFh_calc: logarithm of the constant in front of F^^^ in Eq. (l29i) . 

• sin_chi_calc: sin(x) calculated with Eq. (!20|) . 

• exp_I_omega_chi-calc: e^^^ calculated with Eqs. fll8|20p . 

In cwfcomp.H, the class ODE_integration and member functions extrapolation_in-zero, F_r_u, 
integrationJienrici and operator () performing direct integration of the Coulomb Schrodinger 
equation are defined, as well as the class Coulomh-wave Junctions, with which one can calcu- 
late all Coulomb wave functions. All the routines of the class Coulomhjwavejunctions are in 
cwfcomp.cpp . 

• F_dF_init: initialization of the member variables debut, F_debut and dF_debut used for 
direct integration (see Sec. (14. Sp ). 

• asymptotic-series: calculate the asymptotic series in Eq. (fTOj) for Coulomb wave function 
and derivative. 

• continued Jraction J, continuedJractionJi: calculate the continued fractions of respec- 
tively Eq. dn]) and Eq. 

• F_dF_power_series: calculate Fe^rjiz) and F^^^{z) with Eq. ([H]). 

• asymptotic-expansion-F_dF, asymptotic-expansion_H_dH_scaled: calculate ^£^^(2) and ^(2) 
{Hl^{z) and Hf,^{z) scaled, see Eq. (jSlD) from Eqs. ffT0]l26ll3T|) . 

• F_dF_directJntegration, H_dH_direct_integration\ calculate F^^r^^z) and F^.^{z) {H'^^^{z) 
and Hf^{z)) by direct integration of Eq. ([T]). 

• partiaLderivatives, first_order_expansions: calculate F^^^{z) and F^^^{z) or G(,^^{z) and 
G'^^^{z) with the method of Sec. ^B)- 

• H_dH_from_first_order_expansions: calculate if^^ (2;) and if^^(z) with the method of Sec. ( 12. Sp . 

• H_dH_with_F_dF_and-CF: calculate H'^^^{z) and Hf^{z) with Eqs. (I13|28|29p . 

• H_dH_with_expansion: calculate iJ^^(z) and H'^'^{z) with Eq. (ITTp . 

• F_dF_with-symmetry_relations: calculate F^^y^{z) and F^^{z) for 5R(z) < with Eq. fl2Tl) . 
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Except F_dF_init, all the latter routines are private in the class Coulomhjwavejunctions and 
should not be used directly. 

The following routines of the class Coulomb_wave_functions are public and provide the 
requested Coulomb wave functions: 

• F_dF: calculate Fe^r^i^z) and F^^^{z). 

• G-dG: calculate Gi^j^{z) and G'(^,^{z). 

• H_dH: calculate H^^^tiz) and H^^{z). 

• H_dH_scaled: calculate H^^^{z) and H'^'^{z) scaled (see Eq. ( l3Ti) ). 

The calculation of Gg^r]{,z) and G'(^^^{z) is performed by calculating F^^riiz), Hf^{z), and their 

derivatives, so that one has Gi^rii^) = H^^{z) — iFi^^{z) and G'g^j^{z) = H^^{z) — iF^^^{z). 

The file test_rec_rel.cpp contains additional useful routines using the class Goulomh_wave Junctions: 

• Wronskian_test: function calculating Coulomb wave functions accuracy from Wronskian's 
(see Sec. (1230|) ). 

• cwfJ-tahles-recurrencejrelations: routine calculating Coulomb wave functions with recur- 
rence relations for integer spaced £'s (see Sec. (12.101) ). 

• F_dF_Ltables_rec_reLhelper, cwf_Ltables_rec_reLhelper: routines called by 
cwf-Ltables-recurrence-relations, not intended to be used directly. 

4.2 Use of the program 

Due to its object-oriented programming, cwfcomplex is easy to use. One has to declare first a 
class Goulomh-wavejunctions with three parameters /, eta, and isM-normalized, where / and 
eta are two complex numbers representing i and 1], and isJt_normalized is a boolean equal 
to true if one uses the standard normalization of Coulomb wave functions and false if one 
uses the normalization of Eq. ( l32l) . For example, one declares class Goulomh_wave Junctions 
cwf(isJt-normalized,l,eta);. Then, one can use the member functions of the class cwf. For that, 
one needs the complex value z, the integer omega = ±1 (for H"^^), and two complex numbers 
A, dA to store the Coulomb wave function and its derivative. The instructions to obtain A and 
dA are the following: 

• cwf.F_dF (z,A,dA); to compute F^^r^^z) and F^^^{z). 

• cwf.G_dG (z,A,dA); to compute Gi^rj{z) and G[^^{z). 

• cwf.H_dH (omega,z,A,dA); to compute H^^j^{z) and H'^'^{z). 

• cwf.H_dH_scaled (omega,z,A,dA); to compute H'^^^{z) and H'^'^{z) scaled (see Eq. fl3T])). 
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In order to test the accuracy of previous functions, one has to declare a second class 
Coulomb_wave_functions of parameters isJLnormalized, l+l and eta, for example, 
class Coulomh_wave Junctions cwf_p (isJLnormalized, l+l, eta);. Then, the instruction const 
double W = WronskianJest (z,cwf,cwfjp); provides their relative precision, calculated from 
their Wronskian's and angular momentum recurrence relations (see Sec. f l2.10p ). stored in the 
double W. 

Tables of integer spaced £'s Coulomb wave functions are calculated with the routine 
cwfJdahles-recurrencc-relations. For this, one needs the complex angular momentum of small- 
est modulus l-deh, the number of angular momenta to calculate M, eta, isJt-normalized, the 
number of ^-variables Nz, the one-dimensional complex array of 2;-variables z_tab, six two- 
dimensional Nz X Nl complex arrays F_tab, dF_tab, G-tab, dGdab, Hp_tab, dHp_tab, Hm_tab 
and dHm_tab to respectively store the regular wave functions F and derivatives, the irregular 
wave functions G and derivatives, the outgoing irregular wave functions and derivatives, 
and the incoming irregular wave function H~ and derivatives, so that F_tab[iz][il] will provide 
Fe,n{z) with £ =Ldeb+il and z =zJab[iz] (same for other tables). The instruction is then: 

cwfJJables-recurrencc-relations (Ldeb,Nl,eta,isJt_normalized,Nz,z_tab, 

F_tab, dF_tab, GJ,ab, dGJ,ab,HpJ,ab, dHp_tab,Hm_tab, dHmdab ); 

If one considers a single complex variable z, one can use the following instruction: 

cwfJ-tables_recurrence_relations (Ldeb,Nl,eta,isJt_normalized,z, 

F_tab, dF_tab, GJ,ab, dGJ,ab,HpJ,ab, dHp_tab,Hm_tab, dHm_tab ); 

where F_tab,...,dHm_tab are now one-dimensional arrays of Nl complex numbers {F_tab[il] = 
Fe,n{z) with i =Ldeb+il, same for other tables). 

Examples are provided by the program examples, cpp, which calculates values of different 
Coulomb wave functions and derivatives on a circular z-path of the form z = Re^^, with 
9 E [0 : 27r[, for given Iq, t] and R, their accuracy with the function WronskianJest and a table 
of integer spaced £'s Coulomb wave functions with the routine cwfJJables-recurrencc-relations, 
with the same parameters as before and i starting from Iq. 

4.3 Recommendations 

Even though there are no restrictions for the complex values used in the program, it is advised 
to use 3fJ(£) > — 1. Calculations have indeed been found to be more stable for these values. If 
one has 3?(£) < —1, one can use the symmetry transformation — 1, as Eqs. fll7|18l) 

imply H^^^ = H^_,_^^^ e'^^. 

Due to the direct integration procedures. Coulomb wave functions should be calculated if 
possible for z varying smoothly if one considers tables of Coulomb wave function values. Indeed, 
the complex numbers z, Fi^r^iz), F^^^{z) are stored in the class under the names debut, F_debut 
and dF_debut after each calculation, so that the integration from this point to the next is faster 
and more precise if z varies continuously in the complex plane. Also, it is better for 
to increase on its path as then no continued fraction calculation f{z) (see Eq.(TT^) is needed 
during direct integration. 
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5 Conclusion 



The computation of Coulomb wave functions with all its arguments complex is a very difficult 
task. The single use of power/asymptotic series and continued fractions quickly shows its 
limitation when |?7| or increases. It was found that the range of accessible £, 77, and 

z is greatly augmented by adding the direct integration method of the Coulomb equation. 
Calculations are stable for values of |3(?7)| as important as 80, and can be as large as 

100 as well. This method is particularly stable for the implementation of extremely narrow 
resonant states. However, instabilities appear when both |?7| and are large, and 

smaller or comparable to For example, the values i = 15i, rj = 10, and z = 20e*^ 

used in a calculation similar to the ones presented in Sec. (13. ip provide wrong wave functions 
in the vicinity of 6 = 37r/2. This particular problem can be treated by always accepting the 
value of in F_dF_direct_integration for backward integration (see Sec. (I2.4p ). Other issues 
can be solved by using H_dH_with_expansion in H_dH and H_dH_scaled even if l^l > 1 or 
|Q'(£)| < 1. (The comments in the code beginning with four slashes explain to the user how 
to make modifications accordingly.) Nevertheless, these are solutions for very particular cases 
and cannot be included in a general program. Calculations can also become too long if one 
considers irregular Coulomb wave functions for < |£| < 1 and \z\ < 10~^, as the continued 
fraction of Eq. f|T3|) converges very slowly for z ~ [11] and direct integration cannot be used 
in this region. Even though one encounters numerical problems for large values of |53(£)/3?(£)| 
and |?7| or very small \z\, this program has rendered possible calculations which could not be 
undertaken with standard methods. 
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Test Input 



true 

(1,0.1) 

(50,50) 

100.156 

10 

3 



# Description of the input parameters 
# 

# Boolean: true if ones uses standard normalization, false if one uses alternative normalization. 
7^ Complex: angular momentum 1. 

7^ Complex: Sommerfeld parameter eta. 

# Double: radius R of the path in the complex plane: z = R exp(i theta), theta in [0:2 pi[. 

# Integer: number of points Nz to be considered on the path. 

# Integer: number of points Nl for the recurrence relation : l[rec] = 1+k, k in [0:N1-1]. 

# 

# Compilation: g++ -03 examples. cpp -o run 

# Input instruction: ./run < test.input 

# The output is in test. output . 
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Test Output 



isJt_normalized:true 1:(1,0.1) eta:(50,50) R:100.156 Nz:10 Nl:3 
z:(100.156,0) 

F:(-1.021072923e+15,-2.836755456e+15) F':(1.275057299e+15,-2.729507771e+15) 

G:(2.836755456e+15,-1.021072923e+15) G':(2.729507771e+15,1.275057299e+15) 
H+:(5.673510913c+15,-2.042145845e+15) H+':(5.459015542e+15,2.550114598e+15) 
H-:(7.0774288e-17,1.501204734e-16) H-':(5.671783379e-17,-1.558437769e-16) 
Wronskian test: 1.628444125e-ll 

z:(81.02790609,58.87021973) 

F:(0.01090170509,0.002924757522) F':(0.006665318369,0.003695114571) 
G:(57.24722492,-32.54791917) G':(-42.75162529,10.97359983) 
H+:(57.24430017,-32.53701746) H+':(-42.7553204,10.98026514) 
H-:(57.25014968,-32.55882087) H-':(-42.74793017,10.96693451) 
Wronskian test: 1.307402399e-ll 

z:(30.94990609,95.25401645) 

F:(-2.246133078e-15,2.098754042e-15) F':(-5.747597654e-16,2.506104287e-15) 
G:(-4.367342675e+13,-1.907186698e+14) G':(1.181536987e+14,1.103266317e+14) 

H+:(-4.367342675c+13,-1.907186698e+14) H+':(1.181536987c+14,1.103266317c+14) 
H-:(-4.367342675c+13,-1.907186698e+14) H-':(1.181536987e+14,1.103266317e+14) 
Wronskian test: 2.604836552e-ll 

z:(-30.94990609,95.25401645) 

F:(-3.696304706e-35,8.374503306e-35) F':(5.116568262e-35,9.162544125e-35) 
G:(2.32622983e+33,-4.170545023e+33) G':(2.19801873e+33,4.986604576c+33) 
H+:(2.32622983e+33,-4.170545023e+33) H+':(2.19801873e+33,4.986604576e+33) 
H-:(2.32622983e+33,-4.170545023e+33) H-':(2.19801873e+33,4.986604576e+33) 
Wronskian test: 6.182547672e-12 

z:(-81. 02790609,58.87021973) 

F:(-2.432130956e-67,3.004725207e-66) F':(3.593950134e-66,1.98955822e-66) 
G:(1.065320986e+65,-5.911485321e+64) G':(1.320210869e+64,1.651750319e+65) 
H+:(1.065320986e+65,-5.911485321e+64) H+':(1.320210869e+64,1.651750319e+65) 
H-:(1.065320986e+65,-5.911485321e+64) H-':(1.320210869e+64,1.651750319e+65) 
Wronskian test: 1.818126864e-10 

z:(-100.156,1.22651674e-14) 

F:(7.915510206e-34,-4.070932761e-34) F':(3.191449042e-34,1.291346723e-33) 

G:(4.180561145c+103,8.128671328e+103) G':(-1.326122108e+104,3.277393326e+103) 
H+:(4.180561145c+103,8.128671328e+103) H+':(-1.326122108c+104,3.277393326c+103) 
H-:(4.180561145e+103,8.128671328e+103) H-':(-1.326122108e+104,3.277393326e+103) 
Wronskian test: 1.24249569e-16 
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z:(-81. 02790609,-58.87021973) 

F: (-23318.74764,-17080.04412) F':(28164.502,-34900.64756) 
G:(17080.04412,-23318. 74763) G':(34900.64758,28164.50199) 
H+: (34160.08824,-46637.49527) H+':(69801.29514,56329.004) 
H-:(6.998646405e-06,8.686284079e-06) H-':(1.393963282e-05,-1.022684069e-05) 
Wronskian test: 1.320762168e-ll 

z:(-30.94990609,-95.25401645) 

F: (-3.4196046366+30 ,-3.206140946e+30) F':(4.148673182e+30,-5.870853625e+30) 
G:(3.206140946e+30,-3.419604636e+30) G':(5.870853625e+30,4.148673182e+30) 

H+:(6.412281891e+30,-6.839209271e+30) H+':(1.174170725e+31,8.297346365e+30) 
H-:(4.01631585e-32,5.686566419e-32) H-':(7.772167842e-32,-7.290658234e-32) 
Wronskian test: 1.127990924e-ll 

z:(30.94990609,-95.25401645) 

F:(1.125583254e+40,3.548477279e+39) F':(3.759922307e+38,1.698605313e+40) 
G:(-3.548477279c+39,1.125583254c+40) G':(-1.698605313c+40,3.759922307c+38) 
HH-:(-7.096954559e+39,2.251166509e+40) H+':(-3.397210626e+40,7.519844613e+38) 
H-:(6.373392552e-43,-2.945348841e-41) H-':(-4.036846683e-41,1.270436304e-41) 
Wronskian test: 1.056406572e-ll 

z:(81.02790609,-58.87021973) 

F:(-2.579395538e+32,7.380968215e+32) F':(-9.701448492e+32,1.926734513e+32) 
G:(-7.380968215e+32,-2.579395538e+32) G':(-1.926734513e+32,-9.701448492e+32) 
H+:(-1.476193643e+33,-5.158791075e+32) H+':(-3.853469026e+32,-1.940289698e+33) 
H-:(-4.963907579c-34,-9.779175601e-35) H-':(2.098665903e-34,6.035174252e-34) 
Wronskian test: 1.628966764e-12 



Recurrence relations results for a table of z values. 

z:(100.156,0) 

l[rec]:(l,0.1) 

F:(-1.021072923e+15,-2.836755456e+15) F':(1.275057299e+15,-2.729507771e+15) 
G:(2.836755456e+15,-1.021072923e+15) G':(2.729507771e+15,1.275057299e+15) 
H+:(5.673510913e+15,-2.042145845e+15) H+':(5.459015542e+15,2.550114597e+15) 
H-:(7.0774288e-17,1.501204734e-16) H-':(5.67178338e-17,-1.558437769e-16) 

l[rec]: (2,0.1) 

F:(-9.963131598e+14,-2.756374147e+15) F':(1.235449857e+15,-2.655381893e+15) 
G:(2.756374147e+15,-9.963131598e+14) G':(2.655381893e+15,1.235449857e+15) 
H+:(5.512748294e+15,-l. 992626326+15) H+':(5.310763785e+15,2.470899715e+15) 
H-:(7.256451117e-17,1.54534022e-16) H-':(5.855852326e-17,-1.602324976e-16) 

l[rec]: (3,0.1) 
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F:(-9.584372325e+14,-2.64073598e+15) F':(1.180083758e+15,-2.547127509e+15) 

G:(2.64073598e+15,-9.584372325c+14) G':(2.547127509c+15,1.180083758c+15) 
H+:(5.28147196e+15,-1.916874465c+15) H+':(5.094255019c+15,2.360167515c+15) 
H-:(7.544527529e-17,1.613446827e-16) H-':(6.131348038e-17,-1.670888938e-16) 

z:(81.02790609,58.87021973) 

l[rec]:(l,0.1) 

F:(0.01090170509,0.002924757522) F':(0.00666531837,0.00369511457) 
G:(57.24722493,-32.54791916) G':(-42.75162529,10.97359982) 
H+:(57.24430017,-32.53701746) H+':(-42.7553204,10.98026514) 
H-:(57.25014969,-32.55882087) H-':(-42.74793018,10.9669345) 

l[rec]:(2,0.1) 

F:(0.01069308234,0.002972399334) F':(0.006522115638,0.003689236481) 
G:(57.93801038,-33.60221658) G':(-43.37569649,11.55133447) 
H+: (57.93503798,-33.5915235) H+':(-43.37938573,11.55785658) 
H-:(57.94098278,-33.61290967) H-':(-43.37200725,11.54481235) 

l[rec]: (3,0.1) 

F:(0.01038725067,0.003043902107) F':(0.006311911265,0.003681809291) 
G:(58.95045958,-35.24443848) G':(-44.3055982, 12.46205462) 

H+:(58.94741568,-35.23405123) H+':(-44.30928001, 12.46836653) 
H-: (58.95350348,-35.25482573) H-':(-44.30191639,12.45574271) 

z:(30.94990609,95.25401645) 
l[rec]:(l,0.1) 

F:(-2.246133078e-15,2.098754042e-15) F':(-5.747597651e-16,2.506104288c-15) 
G:(-4.367342673c+13,-1.907186698e+14) G':(1.181536987c+14,1.103266317c+14) 
H+:(-4.367342673e+13,-1.907186698e+14) H+':(1.181536987e-t-14,1.103266317e+14) 
H-:(-4.367342673e+13,-1.907186698e+14) H-':(1.181536987e+14,1.103266317e+14) 

l[rec]:(2,0.1) 

F:(-2.280830781e-15,2.034967775e-15) F':(-6.280365023e-16,2.478500727e-15) 
G:(-4.826447957e+13,-1.907399944e+14) G':(1.213466147e+14,1.081922962e+14) 
H+:(-4.826447957e+13,-1.907399944e+14) H+':(1.213466147e+14,1.081922962e+14) 
H-:(-4.826447957e+13,-1.907399944e+14) H-':(1.213466147e+14,1.081922962e+14) 

l[rec]:(3,0.1) 

F:(-2.331318569e-15,1.939356637e-15) F':(-7.06875848e-16,2.436443775e-15) 
G:(-5.518815356e+13,-1.904610411e+14) G':(1.260154402e+14,1.047559123e+14) 
H+:(-5.518815356e+13,-1.904610411e+14) H+':(1.260154402e+14,1.047559123e+14) 
H-:(-5.518815356e+13,-1.904610411e+14) H-':(1.260154402e+14,1.047559123e+14) 

z:(-30.94990609,95.25401645) 
l[rec]:(l,0.1) 



21 



F:(-3.696304706e-35,8.374503306e-35) F':(5.116568262e-35,9.162544125e-35) 

G:(2.32622983e+33,-4.170545023c+33) G':(2.19801873e+33,4.986604576e+33) 
H+:(2.32622983e+33,-4.170545023e+33) H+':(2.19801873e+33,4.986604576e+33) 
H-:(2.32622983e+33,-4.170545023e+33) H-':(2.19801873e+33,4.986604576e+33) 

l[rec]:(2,0.1) 

F:(-3.976658859e-35,8.274528706e-35) F':(4.8327284e-35,9.349101023e-35) 
G:(2.184541884e+33,-4.231185626e+33) G':(2.351297285e+33,4.898721703e+33) 
H+:(2.184541884e+33,-4.231185626e+33) H+':(2.351297285e+33,4.898721703e+33) 
H-:(2.184541884e+33,-4.231185626e+33) H-':(2.351297285e+33,4.898721703e+33) 

l[rec]:(3,0.1) 

F:(-4.394887926e-35,8.113966164e-35) F':(4.398972708e-35,9.619612834e-35) 
G:(1.968199386e+33,-4.309629699e+33) G':(2.572345582e+33,4.754740671e+33) 
H+:(1.968199386e+33,-4.309629699e+33) H+':(2.572345582e+33,4.754740671e+33) 
H-:(1.968199386e+33,-4.309629699e+33) H-':(2.572345582e+33,4.754740671e+33) 

z:(-81. 02790609,58.87021973) 
l[rec]:(l,0.1) 

F:(-2.432130956e-67,3.004725207e-66) F':(3.593950134e-66,1.98955822e-66) 
G:(1.065320986e+65,-5.911485322e+64) G':(1.320210869e+64,1.651750319e+65) 

H+:(1.065320986c+65,-5.911485322c+64) H+':(1.320210869e+64,1.651750319c+65) 
H-:(1.065320986e+65,-5.911485322e+64) H-':(1.320210869e+64,1.651750319e+65) 

l[rec]: (2,0.1) 

F:(-3.576416599e-67,3.029259022e-66) F':(3.559897692e-66,2.145032835e-66) 

G:(1.030789257c+65,-6.225220366c+64) G':(1.903651657c+64,1.626497749c+65) 
H+:(1.030789257c+65,-6.225220366c+64) H+':(1.903651657c+64,1.626497749e+65) 
H-:(1.030789257e+65,-6.225220366e+64) H-':(1.903651657e+64,1.626497749e+65) 

l[rec]:(3,0.1) 

F:(-5.328192626e-67,3.06145691c-66) F':(3.501151213c-66,2.380004106e-66) 
G:(9.770243181e+64,-6.655971946e+64) G':(2.74025893e+64,1.583938607c+65) 
H+:(9.770243181e+64,-6.655971946e+64) H+':(2.74025893e+64,1.583938607e+65) 
H-:(9.770243181e+64,-6.655971946e+64) H-':(2.74025893e+64,1.583938607e+65) 

z:(-100.156,1.22651674e-14) 

l[rec]:(l,0.1) 

F:(7.915510206e-34,-4.070932761e-34) F':(3.191449041e-34,1.291346722e-33) 
G:(4.180561145e+103,8.128671329e+103) G':(-1.326122108e+104,3.277393326e+103) 
H+:(4.180561145e+103,8.128671329e+103) H+':(-1.326122108e+104,3.277393326e+103) 
H-:(4.180561145e+103,8.128671329e+103) H-':(-1.326122108e+104,3.277393326e+103) 

l[rec]: (2,0.1) 

F:(7.60252401e-34,-4.273651875e-34) F':(3.593368917^34,1.252731129e-33) 
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G:(4.388739394e+103,7.807256555e+103) G':(-1.286466613e+104,3.690136722e+103) 

H+:(4.388739394e+103,7.807256555c+103) H+':(-1.286466613e+104,3.690136722e+103) 
H-:(4.388739394e+103,7.807256555e+103) H-':(-1.286466613e+104,3.690136722e+103) 

l[rec]: (3,0.1) 

F:(7.127789352^34,-4.540648434e-34) F':(4.144057432e-34,1. 192784614^33) 

G:(4.66292605e+103,7.319737507e+103) G':(-1.224905766e+104,4.255655031e+103) 
H+:(4.66292605e+103,7.319737507e+103) H+':(-1.224905766e+104,4.255655031e+103) 
H-:(4.66292605e+103,7.319737507e+103) H-':(-1.224905766e+104,4.255655031e+103) 

z:(-81. 02790609,-58.87021973) 

l[rec]:(l,0.1) 

F:(-23318.74764,-17080.04412) F':(28164.502,-34900.64756) 
G:(17080.04412,-23318. 74763) G': (34900.64757,28164.50199) 
H+:(34160.08824,-46637.49528) H+':(69801.29514,56329.004) 
H-:(6.998646405e-06,8.686284079e-06) H-':(1.393963282e-05,-1.022684069e-05) 

l[rec]:(2,0.1) 

F: (-23227. 13327,-15830.94943) F':(26224.1359,-34846.45541) 
G:(15830.94944,-23227.13327) G':(34846.45543,26224.13589) 
H+:(31661.89887,-46454.26654) H+':(69692.91084,52448.27179) 
H-:(6.890524871e-06,9.17086489e-06) H-':(1.468279939e-05,-1.002406502e-05) 

l[rec]: (3,0.1) 

F:(-22989.33501,-14034.3731) F': (23425.9144,-34604.63695) 
G:(14034.3731,-22989.335) G':(34604.63697,23425.91439) 
H-t-: (28068.7462,-45978.67001) H+':(69209.27392,46851.8288) 
H-:(6.703413121e-06,9.918798258e-06) H-':(1.582831668e-05,-9.679629412e-06) 

z:(-30.94990609,-95.25401645) 
l[rec]:(l,0.1) 

F:(-3.419604636e+30,-3.206140946e+30) F':(4.148673182c+30,-5.870853625e+30) 
G:(3.206140946e+30,-3.419604636e+30) G':(5.870853625c+30,4.148673182e+30) 
H+:(6.412281891e+30,-6.839209271e+30) H+':(1.174170725e+31,8.297346364e+30) 
H-:(4.01631585e-32,5.68656642e-32) H-':(7.772167843e-32,-7.290658234e-32) 

l[rec]:(2,0.1) 

F:(-3.392855448e+30,-3.007724359c+30) F':(3.853513229e+30,-5.788374417c+30) 
G:(3.007724359e+30,-3.392855448e+30) G':(5.788374417e+30,3.853513229e-t-30) 
H+:(6.015448719e+30,-6.785710896e+30) H+':(1.157674883e+31,7.707026457e+30) 
H-:(3.986897101e-32,5.991979415e-32) H-':(8.242309638e-32,-7.310395217e-32) 

l[rec]:(3,0.1) 

F:(-3.339981893e+30,-2.724844995e+30) F':(3.435784336e+30,-5.648469464e+30) 
G:(2.724844995e+30,-3.339981893e+30) G':(5.648469464e+30,3.435784336e+30) 
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H+:(5.449689989e+30,-6.679963785e+30) H+':(1.129693893e+31,6.871568672e+30) 
H-:(3.932323694e-32,6.468436407e-32) H-':(8.977628494e-32,-7.327940398e-32) 

z:(30.94990609,-95.25401645) 
l[rec]:(l,0.1) 

F:(1.125583254e+40,3.548477279e+39) F':(3.759922313e+38,1.698605313e+40) 

G:(-3.548477279e+39,1.125583254e+40) G':(-1.698605313c+40,3.759922313e+38) 
H+:(-7.096954559e+39,2.251166509e+40) H+':(-3.397210626e+40,7.519844627e+38) 
H-:(6.373392558e-43,-2.945348841e-41) H-':(-4.036846684e-41,1.270436304e-41) 

l[rec]:(2,0.1) 

F:(1.092775914e+40,3.208281133e+39) F':(6.877131706e+38,1.638289489e+40) 
G:(-3.208281133e+39,1.092775914e+40) G':(-1.638289489e+40,6.877131706e+38) 
H+:(-6.416562265e+39,2.185551829e+40) H+':(-3.276578978e+40,1.375426341e+39) 
H-:(1.265131507e-42,-3.049938868e-41) H-':(-4.208492895e-41,1.233299224e-41) 

l[rec]: (3,0.1) 

F:(1.043809065e+40,2.740200645e+39) F':(1.099037848c+39,1.550082818c+40) 
G:(-2.740200645e+39,1.043809065e+40) G':(-1.550082818e+40,1.099037848e+39) 
H+:(-5.480401289e+39,2.087618131e+40) H+':(-3.100165636e+40,2.198075696e+39) 
H-:(2.26211583e-42,-3.213069851e-41) H-':(-4.477068961e-41,1.172939493e-41) 

z:(81.02790609,-58.87021973) 
l[rec]:(l,0.1) 

F:(-2.579395538e+32,7.380968215e+32) F':(-9.701448491e+32,1.926734514e+32) 
G:(-7.380968215e+32,-2.579395538e+32) G':(-1.926734514e+32,-9.701448491e+32) 

H+:(-1.476193643e+33,-5.158791075c+32) H+':(-3.853469027c+32,-1.940289698e+33) 
H-:(-4.963907579e-34,-9.779175603e-35) H-':(2.098665903e-34,6.035174253e-34) 

l[rec]: (2,0.1) 

F:(-2.416251979e+32,7.159184131e+32) F':(-9.355630991e+32,1.963184002e+32) 

G:(-7.159184131c+32,-2.416251979c+32) G':(-1.963184002e+32,-9.355630991c+32) 
H+:(-1.431836826c+33,-4.832503957e+32) H+':(-3.926368004e+32,-1.871126198e+33) 
H-:(-5.124854379e-34,-1.067175763e-34) H-':(2.104701743e-34,6.26792112e-34) 

l[rec]: (3,0.1) 

F:(-2.191409619e+32,6.833977857c+32) F':(-8.857660418e+32,2.000904331e+32) 
G:(-6.833977857e+32,-2.191409619c+32) G':(-2.000904331c+32,-8.857660418e+32) 
H+:(-1.366795571e+33,-4.382819238e+32) H+':(-4.001808662e+32,-1.771532084e+33) 
H-:(-5.377076882e-34,-1.205969114e-34) H-':(2.115384163e-34,6.632005233e-34) 



Recurrence relations results for a single z. 

z:(100.156,0) 

l[rec]:(l,0.1) 
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F:(-1.021072923e+15,-2.836755456e+15) F':(1.275057299e+15,-2.729507771e+15) 

G:(2.836755456e+15,-1.021072923c+15) G':(2.729507771c+15,1.275057299c+15) 
H+:(5.673510913e+15,-2.042145845c+15) H+':(5.459015542c+15,2.550114597e+15) 
H-:(7.0774288e-17,1.501204734e-16) H-':(5.67178338e-17,-1.558437769e-16) 

l[rec]:(2,0.1) 

F:(-9.963131598e+14,-2.756374147c+15) F':(1.235449857e+15,-2.655381893e+15) 
G:(2.756374147e+15,-9.963131598e+14) G':(2.655381893c+15,1.235449857e+15) 
H+:(5.512748294e+15,-1.99262632e+15) H+':(5.310763785e+15,2.470899715e+15) 
H-:(7.256451117e-17,1.54534022e-16) H-':(5.855852326e-17,-1.602324976e-16) 

l[rec]:(3,0.1) 

F:(-9.584372325e+14,-2.64073598e+15) F':(1.180083758e+15,-2.547127509e+15) 
G:(2.64073598e+15,-9.584372325e+14) G':(2.547127509e+15,1.180083758e+15) 
H+:(5.28147196e+15,-1.916874465e+15) H+':(5.094255019e+15,2.360167515e+15) 
H-:(7.544527529e-17,1.613446827e-16) H-':(6.131348038e-17,-1.670888938e-16) 
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Table 1: Number of iterations n needed for the continued fraction h^{z) to converge up to 
e = 10"^'^ with Lentz method. The set of parameters is chosen here as £ = 0, 77 = 10, and 
z — X — 2i, with X varying between and 1. 



X 


n 


1 


485 


0.5 


1,679 


0.1 


35,984 


0.05 


137,922 


0.01 


3,146,899 


0.005 


12,107,924 



Table 2: Comparison of the Coulomb wave functions F^^^i^z) and Gg^r^{z) and derivatives cal- 
culated with direct and first-order expansion methods. The second and third columns show 
respectively the real and imaginary parts of the Coulomb wave functions obtained with the 
direct method. The fourth and fifth column provides the relative difference of the same real 
and imaginary parts with those calculated with the first-order expansion respectively. The used 
set of parameters is £ = ilO~^, 77 = 10 -|- ilO~^, and z = 0.1 + ilO~^. 





Direct (real part) 


Direct (im. part) 


Rel. diff. (real part) 


Rel. diff. (im. part) 


F 


4.306-10-^^ 


-9.033-10-1^ 


1.481-10-1" 


6.819-10-1" 


F' 


7.635-10-13 


-1.861-10-1^ 


2.507-10-1° 


5.467-10-1° 


G 


7.787-1011 


1.842-10^ 


3.299-10-1° 


1.053-10-*^ 


G' 


-9.416-1012 


-2.076-10^ 


3.086-10-1° 


2.849-10-1° 
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Table 3: Energies and widths of the 2si/2 proton state of the Hamiltonian h of Eq. ( l36i) as 
a function of the depth of the Woods-Saxon potential Vq given in MeV. 7 denotes the width 
obtained by direct integration of h, and 7c is the width obtained by the approximate current 
formula of Eq. (1371) . Energies are given in MeV and widths in keV. 



Vo (MeV) 


E (MeV) 


7 (keV) 


7c (keV) 


50 


4.510 


6.188-10-1 


6.188-10-1 


51 


3.847 


6.134-10-2 


6.134-10-2 


52 


3.168 


2.597-10-3 


2.597-10-3 


53 


2.477 


2.664-10-5 


2.664-10-5 


54 


1.773 


1.777-10-8 


1.777-10-8 


55 


1.060 


1.260-10-1^ 


1.260-10-1^ 


56 


0.336 


7.738-10-36 


7.738-10-36 


56.4 


4.458-10-2 


4.956-10-121 


4.956-10-121 


56.46 


6.949-10-^ 
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Figure 1: Decimal logarithms of \F£^r]{z)\, \H^^{z)\, and \Hl^{z)\ as a function of the argument 
6 of z, given in units of radian over it. They are respectively represented with full, dashed, and 
dotted lines, z = Rtc^^ with Rt the generalized turning point \rj\ + + 1)| + |?7|^. I and t] 

are here respectively equal to O.li and 50 + O.li. 




I I I 1 I I L_ 

0.5 1 1.5 

(rad % ) 



Figure 2: Same as Fig. ([T]), except that £ and rj are here respectively equal to O.lz and 50 + 501 
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Figure 3: Same as Fig. ([T]), except that i and rj are here respectively equal to O.li and 1 + 80i. 




-100- ,..log,„(IH-(z)l). _ 

I I 1 I I " r" I 

0.5 1 1.5 2 

(rad n ) 

Figure 4: Same as Fig. ([1]), except that £ and r] are here respectively equal to lOOi and 0.1 + O.li 
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Figure 5: Same as Fig. ([1]), except that £ and r] are here respectively equal to 100 + lOi and 
50 + 50i 
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